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Abstract 

Inspired by a recent hyperbolic regularization of Burnett's hydrodynamic equations [A. Bobylev, J. Stat. 
Phys. 124, 371 (2006)], we introduce a method to derive hyperbolic equations of linear hydrodynamics to 
any desired accuracy in Knudsen number. The approach is based on a dynamic invariance principle which 
derives exact constitutive relations for the stress tensor and heat flux, and a transformation which renders 
the exact equations of hydrodynamics hyperbolic and stable. The method is described in detail for a simple 
kinetic model - a thirteen moment Grad system. 

PACS numbers: 05.20.Dd, 51.10.+y 
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I. INTRODUCTION 



Derivation of hydrodynamics from a microscopic description is the classical problem of phys- 
ical kinetics. The Chapman-Enskog method [1] derives the solution from the Boltzmann equation 
in a form of a series in powers of Knudsen number e, where e is a ratio between the mean free path 
of a particle and the scale of variations of hydrodynamic fields. The Chapman-Enskog solution 
leads to a formal expansion of stress tensor and of heat flux vector in balance equations for density, 
momentum, and energy. Retaining the first order term (e) in the latter expansions, we come to the 



Navier-Stokes equations, while next-order corrections are known as the Burnett |2J] (e 2 ) and the 
super-Burnett (e 3 ) corrections flJLO . It has long been conjectured that the inclusion of higher-order 
terms in the constitutive relations for the stress and heat flux may improve the predictive capa- 
bilities of hydrodynamics formulations in the continuum-transition regime where Navier-Stokes 
equations fail. 

However, as it was first demonstrated by Bobylev for Maxwell's molecules [3], even in the 
simplest case (one-dimensional linear deviation from global equilibrium), the Burnett and the 
super-Burnett hydrodynamics violate the basic physics behind the Boltzmann equation. Namely, 
sufficiently short acoustic waves are increasing with time instead of decaying. Bobylev's instabil- 
ity has been also studied by Uribe et al. [4] for hard sphere molecules. This instability contradicts 
the //-theorem, since all near-equilibrium perturbations must decay. This creates difficulties for 
an extension of hydrodynamics, as derived from a microscopic description, into a highly non- 
equilibrium domain where the Navier-Stokes approximation is inapplicable. For example, higher- 
order systems of hydrodynamic equations afforded a better description in certain situations such 
as shock structures on coarse grids, but are prone to small wavelength instabilities when grids 
are refined. Successes and drawbacks of the Burnett computations and a family of extended Bur- 
nett equations aimed at achieving entropy-consistent behavior of the equations have been recently 
summarized in yj]. 

The failure of the CE expansion does not lie in the method itself, but in its truncation to lower 
order levels. This question was studied in some detail for a class of simple kinetic models - 
Grad's moment systems [6] - in Refs. l^flSHIlllll^]. In these works, the Chapman-Enskog 
expansion was summed up exactly which revealed stability of exact hydrodynamics in contrast to 
its finite-order approximations. Alternative ways of approximating the Chapman-Enskog solution 
have been also suggested. 
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Very recently, Bobylev suggested a different viewpoint on the problem of Burnett's hydro- 



dynamics n 1 3M . Namely, violation of hyperbolicity can be seen as a source of instability. We 
remind that Boltzmann's and Grad's equations are hyperbolic and stable due to corresponding H- 
theorems. However, the Burnett hydrodynamics is not hyperbolic which leads to no //-theorem. 
Bobylev [13] suggested to stipulate hyperbolization of Burnett's equations which can also be con- 
sidered as a change of variables. In this way hyperbolically regularized Burnett's equations admit 
the //-theorem (in the linear case, at least) and stability is restored. 

The aim of this paper is to study the issue of hyperbolicity of higher-order hydrodynamics in the 
case where the Chapman-Enskog solution can be found exactly. As a starting point, we consider 
Grad's moment system, linearized at the equilibrium, and assuming unidirectional flow conditions 
(the 1D13M system, according to tw- While simple enough, this model is nontrivial for three 
reasons: 

• Application of the Chapman-Enskog method leads to a rather involved nonlinear recurrent 
relations for the coefficients of the expansion; 

• The Burnett approximation derived from Grad's moment system is identical to the one de- 
rived from the Boltzmann equation for Maxwell's molecules and thus violates hyperbolicity 
and exhibits Bobylev's instability PO; 

HoLllli 



120, and is stable, the question remains whether or not this exact hydrodynamics is 
manifestly hyperbolic. 

The paper is organized as follows: In Sec. HH we derive exact hydrodynamics from the lin- 
earized 1D13M Grad's system. Derivation closely follows [12], and is based on application of 
a dynamic invariance principle which is equivalent to exact summation of the Chapman-Enskog 
expansion. A critical value of the Knudsen number is found beyond which a closed system of 
equations for the locally conserved fields ceases to exist. In Sec. [Ill] we find a class of transforma- 
tions through which exact equations of hydrodynamics can be put in a hyperbolic form, thereby 
answering in affirmative the above question. We also analyze how such transformations affect the 
dissipative nature of the equations. In Sec. [TV] we analyze and compare with conventional and 
earlier approximate solutions provided by (i) the Newton iteration method (Appendix [A]) and (ii) 
Bobylev's hyperbolic regularization of Burnett's equations (Appendix IB]) which turns out to be a 
special case of a more general results presented here. Finally, conclusions are provided in Sec.lVl 



II. HYDRODYNAMICS FROM THE LINEARIZED GRAD SYSTEM 



A. Chapman-Enskog method and Bobylev's instability of Burnett's hydrodynamics 

Point of departure is the linearized Grad's thirteen-moment system in one spatial variable x: 

d t p = -d x u, 

d t u = -d x p - d x T - d x a, 

2 2 

d t T = --d x u - -d x q, (1) 

o t o = --d x u - —o x q a, 

3 15 e 

5 2 

9 t q = ~2 dxT ~ dx ° ~ ^ e q - 

Here p(x, t), u(x, t) and T(x, t) are the reduced deviations of density, average velocity and temper- 
ature from their equilibrium values, respectively, and a(x, t) and q(x, t) are reduced ^-component 
of the nonequilibrium stress tensor and heat flux, respectively. Moreover, e > has a meaning of 
the Knudsen number. The latter is given by the ratio between the mean free path A and a charac- 
teristic dimension of the system L and is the smallness parameter in the Chapman-Enskog method 



yj]. The magnitude of the Knudsen number determines the appropriate gas dynamic regime [14]. 
In a sequel, we use rescaled variables t' = et and x' = ex and omit prime to simplify notation. 

The system (OQ) provides the time evolution equations for a set of hydrodynamic (locally con- 
served) fields [p, u, T] coupled to the nonhydrodynamic fields a and q. The goal is to reduce the 
number of equations in (0Q) and to arrive at a closed system of three equations for the hydrodynamic 
fields only. Thanks to linearity of the system dU) it proves convenient to turn into the reciprocal 
space, and seek for solutions of the form £ = ( k exp(cot + ikx), where £ is a generic function 
p,u,T, a, q, and where A; is a real valued wavenumber. 

Application of the Chapman-Enskog (CE) method to the reduction of the system (OQ) results in 
the following series expansion of the nonhydrodynamic variables into the powers of k: 

oo oo 

^E-rU-E^, (2) 

n=0 n=0 

where the coefficients and qjf 1 are of order ~ k n+1 , q^ ~ k n+1 , and are obtained from 
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a recurrence procedure: 



n-l 



°k 



??° = - <! E ^g^^ + ifor^ |> , (3) 

Lm=0 J 

and where the CE operators d[ m ^ act on the hydrodynamic fields as follows: 

a t Pk 



,{m) I , m = 



m > 1 



-ik(p k + T k ) ,m = 
— ifccr[ m_1 ' ) , m > 1 



d t (m) T fc = ^ 3 k . (4) 

It can be proven that functions a k and have the following structure, for all n = 0, 1, .., 

a k n) = a n (-k 2 ) n iku k ,, 

4 2n+1) = b n (-er +i Pk + c n (-er +i T k , 

1k n) = x n (-k 2 ) n ikp k + y n (-k 2 ) n ikT kl 
Q { k n+1) = z n (-k 2 r +1 u k , (5) 

where a n , . . . , z n are numerical coefficients to be determined. Note the altering structure of ex- 
pansion coefficients of odd and even orders. Substituting © into © and ©, the CE method casts 
into a recurrence equations in terms of the coefficients a n , . . . , z n : 
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o o n n s 

L \ ^ ^ — V O 

O n -|- ~C n ~\- — / C n — m Z m —\ -\- y (X n — m (X m 7r2 n , 

3 3 z — ' 15 

m=l m=0 

n 2 n g 

m=0 m=0 
n „ n q 

EZ \ O 

Q"n—mCm g / ^■n—raVm ^ p. Un+l i 
m=0 m=0 

n 2 " 

m=l m=0 
n ^ n 

Vn+1 %n ^ j g ^ ^ 2/n— mVm C n , 

m=l m=0 

2 2™ " 

Zn+l — %n+l + gl/n+l + ~ ^ ] yn-m z m + ^ ] Z n - m a m — CL n +l- (6) 

m=0 m=0 

System © is solved recurrently subject to the initial conditions, 

4 i 4 2 15 7 

00 = ~3' = ~3' °° = 3' X ° = 2/0 = ~T' Z ° = ~4' 

The initial conditions are obtained by evaluating the functions a k and q k up to the Burnett order 
(see Eq. (fl~5l) below) and identifying the coefficients ao, xq and from the Navier-Stokes ap- 
proximation and the remaining coefficients b , c and z from the Burnett correction. Equation © 
defines six functions, 

oo oo 

A(k) = J2 *n(-k 2 ) n , Zik) = zn(~k 2 ) n . (8) 

n=0 n=0 

Thus, the CE solution amounts to finding functions A, . . . , Z ©. Note that by the nature of the 
CE recurrence procedure, functions A, . . . , Z ([8]) are real- valued functions. Knowing A, . . . , Z 
d8]), we can express the nonequilibrium stress tensor and heat flux as 

a k = ikA(k)u k - k 2 B(k)p k - k 2 C(k)T k , (9) 
q k = ikX(k)p k + ikY(k)T k - k 2 Z(k)u k . (10) 

Upon substituting these expressions into the Fourier-transformed balance equations (0Q), we obtain 
the closed system of hydrodynamic equations which is conveniently written in a vector form, 

a ( x = Mx, (ii) 
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where x = (p k , u k , T^), and the matrix M has the form 

/ 

M = 







\ 



(12) 



-ik 
ik{l-k 2 B) k 2 A -ik(l-k 2 C) 

y * k 2 x -\%k{\-k 2 z) \k 2 v j 

With this, we find the dispersion relation for the hydrodynamic modes u(k) by solving the char 
acteristic equation, 



det (M - ujI) = 0, 



(13) 



with I the unit matrix. 

The standard application of the CE procedure is to approximate functions A, . . . , Z by polyno- 
mials with coefficients found from the recurrence procedure ©. The first non-vanishing contribu- 
tion is the Newton-Fourier constitutive relations, 



a 



(0) 



--iku k ,q^ ] = — -^ikT k . 
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(14) 



which leads to the Navier-Stokes-Fourier hydrodynamic equations. Computing the coefficients 
and q^, we arrive at the Burnett level: 



4 4 2 

Ok = ~-iku k + -k 2 p k - -k 2 T k , 

q k = — —ikT k + -k u k . 



(15) 



The Burnett approximation (fl~5l) coincides with that obtained by Bobylev [3] from the Boltzmann 
equation for Maxwell's molecules. Unlike the Navier-Stokes-Fourier approximation, the Burnett 
constitutive relations (fl~5l) show instability of the acoustic mode, see Fig. [Q 

Thus, the difficulty of the CE method consists in the way the functions A, . . . , Z are approx- 
imated, the standard polynomial approximations lead to unstable hydrodynamic equations. We 
shall now derive closed-form equations for these functions which amounts to summing up the CE 
series exactly. 



B. Invariance Equations 



Summation of the CE series for the functions A, . . . , Z can be done directly from the recur- 
rence relations © following the lines of Ref. [8]. Alternatively but completely equivalently, one 
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FIG. 1: Dispersion relation. Acoustic mode Re(w ac ) for Navier-Stokes and Burnett hydrodynamics. 

n 

can make use of the dynamic invariance principle (DIP) [7]. Here, the set of non-hydrodynamic 
moments {a, q] is still thought in the form © and (flOl) . but the method makes no assumption 
about the power-series representation of the functions A, . . . , Z. The time derivative of {a, q} can 
be computed in two different ways. On the one hand, substituting © and (flOl) into the moment 
system (OQ), we find 

4 8 

d t a = -~iku k - —ikq(X,Y,Z,k) - a(A,B,C,k), 
6 15 

d t q = ~ikT k -ika(A,B,C,k)-^q(X,Y,Z,k). 

On the other hand, the time derivative of {a, q} can be computed due to the closed hydrodynamic 
equations by chain rule: 



d t cr = ^- d t u k + ^- d t p k + ^ d t T k , (16a) 
au k dp k c)T k 

dq dq dq 

d t q = d t u k + ■— d t p k + T^r d tT k - (16b) 
ou k dpk oT k 

Here, the derivatives d t Uk and d t T k are evaluated self-consistently using the functions © and (flOl) 
in the right hand side of The DIP states that the two time derivatives coincide, since the set 
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{a, q} has to solve both the full Grad system and the reduced system. This requirement implies 
a closed set of equations, here referred as invariance equations (IE), relating the six functions 

A(k),..,Z(k): 

— - — A — k 2 (A 2 +B- — +— ) + -k 4 CZ = 0, 
3 V 15 3 ' 3 

—X + B-A + k 2 AB + -k 2 CX = 0, 
15 3 

— Y + C - A + k 2 AC + -k 2 CY = 0, 
15 3 

A + -Z + k 2 ZA -X--Y + -k 2 YZ = 0, 
3 3 3 

k 2 B --X- k 2 Z + k*ZB - -YX = 0, 
3 3' 

----Y + k 2 {C -Z) + k 4 ZC - -k 2 Y 2 = 0. (17) 
2 3 3 

The same equations can be derived upon summation of the CE expansion. Equations (fT7l) are a 

convenient starting point for evaluation of exact hydrodynamics. For k = one recovers the initial 

conditions ©. 



C. Exact hydrodynamic solutions 

The dispersion relation u(k) was found by simultaneously solving numerically the invariance 
equations (fT7l) and the characteristic equation (TT3T) . The resulting hydrodynamic spectrum consist 
of two modes, the acoustic mode, uj ac (k), represented by two complex-conjugated roots of (fT3l) . 
and the real- valued diffusive heat mode, u}&g(k). cf. Fig. [2} 

Among the many sets of solutions {A(k), Z(k)} to the system (TT71 ), the relevant ones are 
continuous functions with the asymptotics: lim^o ^hydr = 0. Remarkably, we find that the solu- 
tion with this asymptotics is unique, and represented by a pair of complex conjugated sets, [S, S*}, 
shown in Figs. [3] and HI Note that a qualitative change of dynamics arises when the diffusive mode 
couples with one of the two non-hydrodynamical modes of Grad's system at a critical wave num- 
ber k c w 0.3023, which is the value where also the Newton method diverges, cf. App. (fAl). By 
the CE perspective, the hydrodynamics of the diffusive mode stops at k c , since, after that point, 
it becomes a complex- valued function coupled with the conjugated non-hydrodynamic mode, see 
Fig.|2l Essentially, for k > k c , the CE method does not recognize any longer the resulting diffusive 
branch as an extension of a hydrodynamic branch. Also, the set of solutions [S, S*], real valued for 
k < k c , continues upon a complex manifold, cf. Fig. [3] We notice that the occurrence of a pair of 
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FIG. 2: (Color online) Dispersion relation for the linearized 1D13M Grad system (Q]). The unique solution 
of hydrodynamical modes obtained from ( fT3l ) with ( fTTl ) coincides with the real parts of the modes of the 
original Grad system (the plot also shows when pairs of conjugate complex roots appear), the solution of 
the original system £T|) features five u's, while the exact solution of (TT3T ) with (fTTT ) has three lo's for each k 
and degenerated over the hydrodynamic branches at k > k c . 

complex conjugated sets of solution is very plausible due to symmetry reasons: inserting S into the 
dispersion relation, we obtain a pair of complex conjugated acoustic modes [cu^S, k),ul c (S, k)] 
plus one of the complex modes resulting from the extension of the diffusive branch for k > k c ; 
whereas, through S*, we obtain, symmetrically, the two latter conjugated modes, plus one of the 
conjugated acoustic modes. 

As a further evidence of this close coupling, we also notice the occurrence of an intersection 
between the real parts of the hydrodynamical modes Re(o; ac ) and Re(a>diff) after the critical point, 
at k = k coup \ 0.383. Therefore the message extracted from the study of Grad's system (OQ) is 
that the set of hydrodynamic equations for [p, u, T] provides, as expected, stable solutions, when 
taking into account all the orders of CE expansion - which corresponds to solving the system of 
invariance equations (fTTT) . And, there is no closed set of hydrodynamic equations after k c , even 
though the acoustic mode extends smoothly beyond k c , as is visible in Fig. [2l 

Thus, the exact hydrodynamics as derived by the summation of the CE expansion (or, equiv- 
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FIG. 3: (Color online) Imaginary parts of coefficients AtoZ solving (fTTT ). Shown is the unique solution 
leading to hydrodynamic branches, cf. Fig. [2l which is symmetric about the real axis. 

alently, from the invariance equations) extends up to a finite critical value k c . No stability viola- 
tion occurs, unlike in the finite-order truncations thereof. While we have evaluated the functions 
A, . . . , Z numerically, two questions remained open: 

• Is the (stable) exact hydrodynamics also hyperbolic? 

• If so, how to retain hyperbolicity in the approximations? 
In the next section we shall answer the first of these questions. 

HI. HYPERBOLIC TRANSFORMATION FOR EXACT HYDRODYNAMICS 

Equation (Q]) for the Fourier component vector x = (p k , u k , T k ) reads d t x = Mx with M from 
(TP2l) . By explicitly re-introducing the Knudsen number e, i.e. by replacing k by ke in M and 
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FIG. 4: (Color online) Real parts of complex-valued functions A,...,Z solving (fT71) . It is clearly visible 
that the solution matches the initial condition (0). 

further distinguishing between the real and imaginary matrix elements in M, we can write 

<9 t x = [Re(M) — iIm(M)]x, (18) 

oo 

Re(M) = ^(-l) n R (n) £ 2n+1 = eR (0) -e 3 R (1) + O(£ 5 ), 

n=0 

oo 

Im(M) = ^(-l)«l(™)e 2 ™ = - £ 2 lW + £ 4 I( 2 ) - 0(e 6 ), (19) 

n=0 

rearranged such that the Knudsen number expansion coefficients become visible. We find that 
the operators Re(M) (real part) and Im(M) (imaginary part) involve the following real-valued 
operators (for all n > , i.e., with the convention a_i = c_i = Z-\ = 1 and Kronecker symbol 
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si 



j(n) _ k 2n+l 



( 



p/n) = £,2n+2 



5n,0 

6 n _i c n _i 

§2 n -l 



a n 

2„ n 2, 



(20) 



y gX n 3?/n / 

Equations of hydrodynamics (TT8T) are hyperbolic and stable provided that we can find a trans- 
formation of hydrodynamic fields such that 

(i) Re(M) and Im(M) are both real and symmetric, and 

(ii) Re(M) has negative semidefinite eigenvalues. 

Therefore, we seek a transformation z = Tx which produces a symmetric matrix M' = TMT 1 
and we wish to see if Re(M') = Re(TMT _1 ) is negative semidefinite. We consider the equations 
of exact hydrodynamics, i.e. equations (fT8l) provided that functions A, .., Z (|2Q|) are solutions to 
the invariance equations (TT71) . After a few algebra which we do not recapitulate here, we obtain a 
particular transformation matrix T which solves the problem. It is a member of a whole class of 
effectively equivalent transformations, and can be written as 

( T pp T p t ^ 

_ 1 



V 



T 





T p t 





T 














(21) 



/ 



with the nonvanishing components 



- pp 



y /3X + 2Y[[Z]y 

T uu = y/X[[3B - 2Z[[C}} - 2C]] + 2Y[[B]}[[Z}1 
3[[C)]X 



>3X + 2Y[[Z]] 



y/3[[C\] (Y[[B}) - [[C))X), 
where we have introduced the following symbolic notation: 

[[•]] = 1 " (kef . . 
13 



(22) 



(23) 



The transformation T (|2TI) symmetrizes M and renders the system hyperbolic, as can be verified 
by straightforward computation of M' from (1271) . (I25bl) . We further notice, that T contains only 
even powers of (ke) because the same is true for the coefficients A-Z. 

Next we calculate the eigenvalues Ai 2 ,3 of Re(M') - containing transport coefficients - to 
obtain a remarkably simple result: 

Ai =0, A 2 = k 2 eA, A 3 = \k 2 eY. (24) 



From the analysis of the previous section, it follows that the nontrivial eigenvalues (1241) are nega- 
tive semidefinite for all ke (see Fig. |4] which displays the exact numerical solutions for A and Y). 
Hence, the equation describing hyperbolic hydrodynamics (also hyperbolic up to an arbitrarily 
selected order e n , a feature to be used in the next section) attains the form: 

d t z = M'z, (25a) 
M' = TMT 1 (25b) 

for the vector z = {p k , u k , T k } = Tx of transformed hydrodynamic variables, and where M' is 
symmetric and has seminegative eigenvalues. To summarize, 

Hyperbolicity: (M') T = M', (26a) 
ty ■ f -t / Tr t Re ( M ')] ^ °> 

Dissipativity: < (26b) 
[ det[Re(M')] > 0. 

Equation (|25l) with (|2TI) and (fT2l) satisfying (|26l) is the main result of this paper. The occurrence 
of negative eigenvalues in the exact solutions, together with the existence of a transformation T 
which makes the equations hyperbolic, proves that exact hydrodynamics (OQ), without approxima- 
tions, is stable. In the remainder of this paper we shall make use of the hyperbolicity of exact 
hydrodynamics in order to establish approximate hydrodynamic equations which retain this prop- 
erty. 

IV. LOWER ORDER HYPERBOLIC AND STABLE HYDRODYNAMICS 
A. Approximations on the hyperbolic equations 

In applications, one is interested in using truncated hydrodynamic equations by taking into 
account only lower orders of the Knudsen number e. In this case, the functions A,...,Z are 
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replaced by their lower-order approximations, and they can be generally written - as shown already 
in Eq. © - as polynomials truncated to an arbitrary order n. Their coefficients are usually derived 
through the CE recurrence equations, as outlined above. With the exact numerical solution at hand, 
we can also find, at any given order of approximation, the optimal interpolating functions A, .., Z 
solving (TT71) , a method we wish to recommend, and which has been worked out in Tab. HI Exact 
hydrodynamics, as described by Grad's system (OQ), terminates at k c . In this regime we can perform 
a Taylor expansion, up to any order n, upon the elements of all the three matrices T, M, and T _1 . 
Thus, the approximations are done on the manifestly hyperbolic equation (1251) in such a way as 
to retain hyperbolicity and stability in each order of approximation. Is is worthwhile noticing that 
the eigenvalues, upon approximating Eq. (1251) to a polynomial order n, transform in a canonical 
manner: 



and, depending upon the polynomial coefficients, and in particular depending on the sign of the 
highest order coefficients a n , y n , the eigenvalues A2,3 diverge to ±00 for ke — > 00, but stay negative 
for k < k c , if we use coefficients according to the method summarized in Tab. HI We shall now 
consider a few examples of the suggested procedure. 

B. Euler and Navier-Stokes equations 

For the zeroth-order term, Im(M) = I^ ^ (Euler), the transformation is, according to (1221) . 
given by a diagonal matrix with entries T pp = T uu = 1 and Ttt = a/3/2, all eigenvalues are 
identically zero. The first order term, linear in the Knudsen number (Navier Stokes) is obviously 
stable as well; all eigenvalues are negative semidefinite since a = —4/3 and y = —15/4 are both 
negative. 

C. Hyperbolic regularization for the Burnett level 

The Burnett equations are unstable without regularization. For this level of description, second 
order in the Knudsen number e, with Im(M) = — e 2 !^ 1 ', upon inserting the required exact so- 
lutions at vanishing wave number, a , .., z from (fTTT) . cf. Tab. HI into (121]) . ([22]) . the transformation 




(27) 



m=l 



m=l 
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matrix achieving a symmetric Im(M') reads 



/ 



1 + l{ke) 



2 







\ 



T 







1 







(28) 



V 












29 
8^ 



{kef j 



This transformation coincides with the one derived by Bobylev's hyperbolic regularization 
method fl^l . specified for the present model (an alternate derivation which follows closely Ref. 
J13I is given in Appendix |B]). Notice, that up to the Burnett level only the polynomial coefficients 
at vanishing wave number, listed in the first row of Tab. U enter the transformation T, which can 
be indirectly also inferred from the eigenvalues, cf. Eq. (|27T) . 

D. Beyond the Burnett level 

In Tab. HI we provide not only coefficients, but also ranges of applicability for the given co- 
efficients of A — Z which can be used if we extend the procedure to higher order. The optimal 
coefficients are provided by the least squares fit of the numerical data for exact hydrodynamics, 
see Tab. H Within the given ranges, the eigenvalues of Re(M') are negative semidefinite, i.e., the 
spectrum of the acoustic mode to^k) of the corresponding hyperbolic hydrodynamic system is 
then stable for all wavelengths. 

E. Application: Hyperbolic regularization for the super-Burnett level 

Finally, in order to present explicit illustration of the approximation strategy presented in 
Sec. Hill we present the equations on the next, Super-Burnett, level, which takes into account 
terms up to the order (ke) 3 . The equations of change for the transformed variables z read 



d t z 



M'z 



(29) 
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method 


apply at n 


a n /n\ 


b n /n\ 


c n /n\ 


x n /n\ 


Vn/nl 


z n /n\ 





ke < 0.03 


-4/3 


-4/3 


2/3 





-15/4 


-7/4 


1 


ke < 0.17 


-4/3 


-4/3 


2/3 





-15/4 


-7/4 




1 


1.132 


2.536 


-3.735 


-0.716 


5.873 


9.953 


2 


ke < 0.25 


-4/3 


-4/3 


2/3 





-15/4 


-7/4 




1 


0.706 


1.156 


-2.500 


0.309 


4.652 


7.053 




2 


-1.304 


-4.095 


3.720 


3.030 


-3.741 


-8.902 


3 


ke < 0.28 


-4/3 


-4/3 


2/3 





-15/4 


-7/4 




1 


1.104 


3.042 


-4.055 


-1.123 


5.903 


9.718 




2 


0.329 


3.669 


-2.678 


-2.861 


1.398 


2.023 




3 


0.648 


3.083 


-2.540 


-2.340 


2.040 


4.333 



TABLE I: Polynomial coefficients introduced in ([8]) obtained from the exact numerical solution, cf. Fig. [4] 
by requiring that deviations between exact and polynomial series at given order of the method (first column) 
stay below 1% (i.e. would be invisible in the plot). We used the symmetrized functions A(k) + A(— k) over 
the whole real axes for k when performing the fits in order to enforce correct symmetry. This criterion 
corresponds to a regularization procedure which produces stable results up to the limit (ke) < (ke) c = 
0.3023, as is easily verified, and leads to a recommended range of (high precision) applicability of the 
method (second column). For convenience we list faculty-rescaled series coefficients. These coefficients 
are essentially the prefactors for higher order correction terms in hydrodynamic equations and can be used 
to study the intermediate Knudsen number regime < fee < {ke) c . As described in the text part, with a 
suitable transformation matrix T these choices lead to very convenient hyperbolic differential equations for 
the hydrodynamic fields x = (p, u, T). 
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and transformation 
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where third order terms are not present because T is symmetric in k. The Burnett level (|28l) . where 
X\ disappears, is immediately recovered from (f3T1) . The hydrodynamic variables are restored using 
x = T _1 z, with the inverse transformation (suitable at the super-Burnett level), which, due to our 
(arbitrarly) chosen normalization factor T uu in (|2T|) only slightly differs from T: 
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To complete the "simulation algorithm" using (129b — (1321). we need numerical values for x\ and y\, 
and an initial condition for x, or z. One solves the hyperbolically stable system for z, and finally 
calculate x via T -1 . Suitable values for the coefficients are those given in Tab. H for method 1: 
yi = 5.873 and x 1 = —0.716, because higher order coefficients such as x 2 do not enter. The 
equations of this section should allow to study the regime < ke < 0.17 very accurately. For the 
remaining regime, 0.17 < ke < (ke) c the presented equations are also stable and hyperbolic, but 
not as accurate. They are, by definition, more accurate compared with the ones obtained using the 
recursion method. The equations offered in this section serve as an example on how to use our 
more general result, Eq. (|2T|) . 



V. CONCLUSIONS 



In this paper, we have considered derivation of hydrodynamics for a simple model (OQ) for which 
- as we have demonstrated - all details can be explicitly studied. The main finding is that the exact 
hydrodynamic equations (summation of the Chapman-Enskog expansion to all orders) are mani- 
festly hyperbolic and stable. To the best of our knowledge, this is the first complete answer of the 
kind. We have also suggested a way of approximating the higher-order hydrodynamic equations 
using accurate numerical solution of the invariance equations and expansion of the transformation 
which renders the system hyperbolic. The study supports the recent suggestion of Bobylev on the 
hyperbolic regularization of Burnett's approximation, and reduces to the latter in a special case. 
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We conclude this paper with a few comments on the possible extensions of the present ap- 
proach, (i) The technique of deriving exact hydrodynamics and/or hyperbolic approximations 
thereof can be readily applied to linearized Grad's systems with a larger number of moments. In 
particular, we were able to extend the present derivation to the three-dimensional thirteen moment 
system, the results qualitatively agree with the one-dimensional case considered above and will 
be reported separately, (ii) It is also possible to apply the present techniques to derive exact hy- 
drodynamics from the dynamically corrected Grad's systems, first introduced in fl 1 SM and studied 
in some detail in [11611 . The latter equations have arguably better properties than the Grad's equa- 
tions, especially in the moderate Knudsen number regime where the linearized systems become 
relevant to study of micro-flows, (iii) In this paper, we were addresing the boundary conditions for 
neither the Grad's systems nor for the higher-order hydrodynamic equations. As is well known, 
this question remains essentially open for both. Therefore, a different and intriguing field of ap- 
plications of the present technique is the recently established lattice Boltzmann hierarchy (LBH) 



L17, 



18 



19, 



20, 



21, 



221, 



23, 



24, 



25(1 • Although the primitive variables in the LBH are populations 
of carefully chosen discrete velocities, the LBH equations can be cast into a form of moment 
systems similar to Grad's equations. The crucial advantage of the LBH above Grad's systems is 
that the former is equipped with pertinent boundary conditions derived directly from the Maxwell- 
Boltzmann theory II 19m . Recently, it has been demonstrated, both numerically and analytically, that 



the LBH is capable of capturing such phenomena as slip and nonlinear Knudsen layers [24, 
The present techniques can be applied for reducing higher-order lattice Boltzmann models with 
advantage for the numerics. However, this goes beyond the scope of this paper, interested reader 



is directed to [23 



24 



250 for details. 
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FIG. 5: Dispersion relations uj{k) for acoustic and diffusive modes obtained via Newton iteration. In the 



plots, shown is also the approximation obtained through Bobylev's hyperbolic regularization (HR) [13]. 
Newton iterations fail for k > k c = 0.3023. 

APPENDIX A: NEWTON ITERATION 

The analytical complexity of either the CE method or the invariance equations is overwhelming 
when we regard systems other than the linearized Grad system, such as the Boltzmann equation. 
Approximate solutions are, then, the only feasible approach. In this section we shall describe 
application of the Newton iteration method to the invariance equations. We used Newton's method, 
cf. Fig. [51 to solve iteratively the equations (OQ), taking, as initial condition, the Euler approximation 
(corresponding to a non-dissipative hydrodynamics: A = ... = Z = 0), which leads, after the 
first iteration, to the same result achievable, alternatively, through a technique of partial summation 
I2] of the CE expansion: essentially, a sort of regularized Burnett approximation. It is seen in Fig. [5] 
that Newton iterations converge rapidly to the exact hydrodynamics in the domain of its validity, 
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APPENDIX B: BOBYLEV'S HYPERBOLIC REGULARIZATION 



This appendix reviews a recent approach by Bobylev ft 1 311 and establishes a connection to the 
second-order variant of our approach. We use the original notation of Ref. Il3ll to facilitate com- 
parisons. 

As was demonstrated in Ref. [| 1 3M . after truncating the CE expansion at the Burnett level, the 
(linearized) equation of hydrodynamics takes the general form: d t x + i(B + e 2 B 1 )x + eAx + 
0(e 3 ) = 0, where x is the vector of hydrodynamics variables [p, u, T] and the operators B , A and 
Bi refer, respectively to Euler, Navier-Stokes and Burnett level of approximation. B = B + e 2 Bi 
is a real non-symmetric operator for e > 0. When applied to the Grad's system (Q~|), these findings 
are a special cases of equation (fT9b with (1201) upon identifying B n = (— l) n I^ n \ A = A and 
A n = (— l) n Rw. The loss of symmetry of the operator B, was identified as the reason of the 
instability occurring in the Burnett equations. In order to cure this loss of symmetry, HR introduces 
a symmetric a real valued operator R and defines a change of variables such that z = x + e 2 Rx - or 
in our notation above, T = (1 + e 2 R) . Hence, the resulting equation of hydrodynamics attains the 
form: z t + i[B + e 2 (B l + RB - B R)}z + eAz + O(e 3 ) = 0, more generally z = TMT x z. The 
suggested regularization consists in writing T _1 as a polynomial (Taylor) expansion in powers of 
e and in truncating it, as for T, at second order. 

HR, therefore, provides a regularization which is exact up to the order e 2 . The operator R 
has to be chosen in such a way that B x = B\ + [R, B ] is real and symmetric, where [R, B ] = 
RB — B R. It is instructive to consider the HR as applied to the example of the 1D13M (Q~|) which 
has originally been written, in matrix notation, as: 
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Expression (IB 1 1) offers those first terms of (1191) with (1201) for which the coefficients (a = —4/3, .. , 
z = —7/4) are analytically known, cf. Tab.Hfor all values. Equation (|B1|) can hence equivalently 
be formulated as M = eH^ — — e 2 l^) + 0(e 3 ). To apply the regularization procedure 
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to the system ©, one needs to make matrix B symmetric (it corresponds to restoring the hy- 
perbolicity of Euler equations, through a transformation T a ). Then, introducing a real-valued, 
symmetric (diagonal) matrix R with diagonal elements a(k), b(k), and c(k) (which corresponds 
choosing a diagonal T), and imposing the symmetry of the resulting operator B 1 (more generally, 
of TI^Tr 1 ), the coefficients are interrelated as follows 111311 



? 2Q 
a(k) = b(k) + -k\ c(k) = b(k) - -k\ 



(B2) 



Notice, the transformation R R is a special case of Eq. (|2Tj) . The resulting operator B 1 is given 
by: 



B x = Bt + [R, B Q ] = Br + b(k)[I, B ] + [m ih B ) 
= B x + [iriy.Bo], 



(B3) 



and therefore unique (independent of b(k)). Hence, the hydrodynamic equations resulting from 
HR as applied to 1D13M attain the form: 
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Since Eq. (IB4I) is a special case of the more general Eq. (|30l) . the connection to Bobylev's work 
has been explicitly established. 
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